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Our goal is to carry out high-precision nuclear structure calculations in connection 
with Radioactive Ion Beam Facilities. The main challenge for the theory of drip line 
nuclei is that the outermost nucleons are weakly bound (implying a large spatial 
distribution) and that these states are strongly coupled to the particle continuum. 
For these reasons, the traditional basis expansion methods fail to converge. We 
overcome these problems by representing the nuclear Hamiltonian on a lattice 
utilizing the Galerkin method with Basis-Spline test functions. We discuss tests 
of the numerical method and applications to the deformed shell model, HF+BCS 
and HFB mean field theories. 



1 Introduction 

In recent years, the area QfjJiuclear structure physics has shown substantial 
progress and rapid growth W. With detectors such as GAMMASPHERE and 
EUROGAM, the hmits of total angular momentum and deformation in atomic 
nuclei have been explored, and new neutron rich nuclei have been identified in 
spontaneous fission studies. Gamma-ray detectors under development such as 
GRETA a will have improved resolving power and should allow for the iden- 
tification of weakly populated states never seen before in nuclei. Particularly 
exciting is the proposed construction of a next-generation ISOL FACILITY 
in the United States which has been been identified in the 1996 DOE/NSAC 
Long Range Plan u as the highest priority for major new construction. 

These experimental developments as well as recent advances in computa- 
tional physics have sparked renewed interest in nuclear structure theory. In 
contrast to the well-understood behavior near the valley of stability, there are 
many open questions as we move towards the proton and neutron driplines 
and towards the limits in mass number (superheavy region) . While the proton 
dripline has been explored experimentally up to Z—83, the neutron dripline 
represents mostly "terra incognita". In these exotic regions of the nuclear 
chart, one expects to see several new phenomena: near the neutron dripline, 
the neutron-matter distribution will be very diffuse and of large size giving 
rise to "neutron halos" and "neutrons skins" . One also expects new collec- 
tive modes associated with this neutron skin, e.g. the "scissors" vibrational 
mode or the "pygmy" resonance. In proton-rich nuclei, we have recently seen 
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both spherical and deformed proton emitters; this "proton radioactivity" is 
caused by the tunnehng of weakly bound protons through the Coulomb bar- 
rier. The investigation of the properties of exotic nuclei is also essential for 
our understanding of nucleosynthesis in stars and stellar explosions (rp- and 
r- process). Our primary goal is to carry out high-precision nuclear structure 
calculations in connection with Radioactive Ion Beam Facilities. Some of the 
topics of interest are the effective N-N interaction at large isospin, large pair- 
ing correlations and their density dependence, neutron halos/skins, and proton 
radioactivity. Specifically, we are interested in calculating observables such as 
the total binding energy, charge radii, densities pp,„(r), separation energies 
for neutrons and protons, pairing gaps, and potential energy surfaces. 

There are many theoretical approaches to nuclear structure physics. For 
lack of space, we mention only three of these: in the macroscopic - microscopic 
method, one combines the liquid drop / droplet model with a microscopic 
ell correction frorELa deformed single-particle shell model (MoUer and Nix 
Nazarewicz et al. u) . For relatively light nuclei, it is possible.to diagonalize 
the nuclear Hamiltonian in a shell model basis. Barrett et al. □ have recently 
carried out large-basis no-core shell model calculations for p-shell nuclei. A 
different approach to the interacting nuclear shell model i&^the Shell Model 
Monte Carlo (SMMC) method developed by Dean et al. El which does not 
involve matrix diagonalization but a path integral over auxiliary fields. This 
method has been applied to fp-shell and mediunaJieavy nuclei. Finally, for 
heavier nuclei one utilizes either nonrelativistic BclM or relativistic I11H13 mean 
field theories. 

2 Outline of the theory: HFB formalism in coordinate space 

As we move away from the valley of stability, surprisingly little is known about 
the pairing force: For example, what is its density dependence? Large pair- 
ing correlations are expected near the drip lines which are no longer a small 
residual interaction. Neutron-rich nuclei are expected to be highly superfluid 
due to continuum excitation of neutron "Cooper pairs" . The Hartree-Fock- 
Bogoliubov (HFB) theory unifies the HF mean field theory and the BCS 
pairing theory into a single selfconsistent variational theory. The main chal- 
lenge in the theory of exotic nuclei near the proton or neutron drip line is 
that the outermost nucleons are weakly bound (which implies a very large 
spatial extent), and that the weakly-bound states are strongly coupled to the 
particle continuum. This represents a major problem for mean field theories 
that are based on the traditional shell model basis expansion method in which 
one utilizes bound harmonic oscillator basis wavefunctions. As illustrated in 
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Figure 1. Inadequacy of shell model basis near the drip lines; need for high-accuracy lattice 
representation. 



Figure |l| a weakly bound state can still be reasonably well represented in the 
oscillator basis, but this is no longer true for the continuum states. In fact, 
Nazarewicz et al. □ have shown that near the driplines the harmonic oscilla- 
tor basis expansion does not converge even if iV = 50 oscillator quanta are 
used. This implies that one either has to use a continuum-shell model basis 
or one has to solve the problem directly on a coordinate space lattice. We 
have chosen the latter method. _ 

Several years ago, Umar et al. E3 have developed a three-dimensional 
HF code in Cartesian coordinates using the Basis-Spline discretization tech- 
nique. The program is based on a density dependent effective N-N interac- 
tion (Skyrme force) which also includes the spin-orbit interaction. The code 
has proven efficient and extremely accurate; it incorporates BCS and Lipkin- 
Nogami pairing, and various constraints. The configuration space Hartree- 
Fock approach has had great successes in predicting systematic trends in the 
global properties of nuclei, in particular the mass, radii, and deformations 
across large regions of the periodic table. 
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So far, our attempts to generalize this 3D code to include self-consistent 
pairing forces (Hartree-Fock-Bogoliubov theory on the lattice) have proven too 
ambitious. The reason may be the lack of a suitable damping operator in 3D. 
We have therefore taken a different approach and developed a new Hartree- 
Fock + BCS pairing code in cylindrical coordinates for axially aattmetric 
nuclei, based on the Galerkin method with B-Spline test functions olla. The 
new code has been written in Fortran 90 and makes extensive use of new data 
concepts, dynamic memory allocation and pointer variables. Extending this 
code, we believe that it will be easier to implement HFB in 2D because one can 
use highly efficient LAPACK routines to diagonalize the lattice Hamiltonian 
and does not necessarily rely on a damping operator. 

We outline now our basic theoretical approach for lattice HFB. As is 
customary, we start by expanding the nucleon field operator into a complete 
orthonormal set of s.p. basis states 0i 

^t(r,s)^^ cUUr,s) (1) 

i 

which leads to the Hamiltonian in occupation number representation 

H = < i\ t \j > cl Cj + ^ < ij\ \mn > c\ c] c„ c™ . (2) 

Like in the BCS theory, one performs a canonical transformation to quasipar- 
ticle operators $, $^ 



i3 \ _ f W V'< \ f c 



(3) 



The HFB ground state is defined as the quasiparticle vacuum 

Pk\^o>= . (4) 
The basic building blocks of the theory are the normal density 

p,^ =< ci>o\c] c,\^o >= iV*V^),j (5) 
and the pairing tensor 

K,, =< ^o\cj c,\<Po >^ iV*U^hj (6) 
from which one can form the generalized density matrix 

= ( -f,* I - p*) ^ = 7^2 = 7^ . (7) 
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Using the definition of the HFB ground state energy 

£:'(7^) =< $o|^f - A7V|$o > (8) 
we derive the equations of motion from the variational principle 

d [E'{n) - tr A(7^2 - 7^)] = (9) 
which results in the standard HFB equations 

[H,7^] = (10) 
with the generalized single-particle Hamiltonian 

H - y, dE'/dp, A = dE'/dK* . (11) 

Our goal is to transform to a coordinate space representation and solve 
the resulting differential equations on a lattice. For this purpose, we define 
two types of quasi-particle wavefunctions , (j)2 

</>t(i?„,m) = ^C/„ , (/.2(i?„,m) 0,(ra) (12) 

i i 

in terms of which the particle density matrix for the HFB ground state as- 
sumes a very simple mathematical structure □ 



Po 



(r,a,r',a') = < $o| ^^(rV) i^(m) |$o > 



= 0,(ra) 0*(rV') = ^ MEn,ra) <t>l{En,v' a') . (13) 

In a similar fashion we obtain for the pairing tensor 

KO (r,cr,r',cr') = < $o| ?A(r'cr') ?/>(rcr) |$o > 

OQ 

= Y,^^3 <^'M </'j(r'^^') = MEn,ra) 0t(S„,rV') . (14) 

i,j E„>0 

For certain types of effective interactions (e.g. Skyrme forces) the HFB 
equations in coordinate space are local and have a structure which is reminis- 
cent of the Dirac equation □ 



h -{h-\) J V02(r) ; V'^2(r) ; ' 

where h is the "particle" Hamiltonian and h denotes the "pairing" Hamilto- 
nian. 
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The various terms in h depend on the particle densities Pq(v) for protons 
and neutrons, on the kinetic energy density Tq(r), and on the spin-current ten- 
sor Jij(r). The pairing Hamiltonian h has a similar structure, but depends on 
the pairing densities ^^(r), Tq(r) and Jij{r) instead. Because of the structural 
similarity between the Dirac equation and the HFB equation in coordinate 
space, we encounter here similar computational challenges: for example, the 
spectrum of quasiparticle energies E is unbounded from above and below. 
The spectrum is discrete for \E\ < —A and continuous for \E\ > —A. In the 
case of axially symmetric nuclei, the spinor wavefunctions 01 (r) and (/'2(r) 
have the structure 

3 Computational method: Spline-Galerkin lattice 
representation 

For nuclei near the p/n driplines, we overcome the convergence problems 
of the traditional shell-model expansion method by representing .the nuclear 
Hamiltonian on a lattice utilizing a Basis-Spline expansion tSllalia. B-Splines 
Bf^{x) are piecewise-continuous polynomial functions of order (Af — 1). They 
represent generalizations of finite elements which are B-splines with M = 2. 
A set of fifth-order B-Splines is shown in Figure |[ 

Let us now discuss the Galerkin method with B-Spline test functions. We 
consider an arbitrary (differential) operator equation 

Ofix)-gix)=0. (17) 

Special cases include eigenvalue equations of the HF/HFB type where O = h 
and g{x) = Ef{x). We assume that both f(x) and g(x) are well approximated 
by Spline functions 

fix) « fix) ^ Yl B^i^y ' ~9{^) « 9{^) ^ E ^*'(^)^' • (18) 

i=l i=l 

Because the functions fix) and gix) are approximations to the exact functions 
fix) and gix), the operator equation will in general only be approximately 
fulfilled 

Ofix) - gix) = Rix) . (19) 

The quantity Rix) is called the residual; it is a measure of the accuracy of 
the lattice representation. We multiply the last equation from the left with 
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Figure 2. Set of fifth-order B-Splines for fixed boundary conditions. 



the spline function Bk{x) and integrate over x 



v{x)dxBk{x)0 f {x) - v{x)dxBk{x)g{x) = v{x)dxBk{x)R{x) . (20) 



We have included a volume element weight function v{x) in the integrals to 
emphasize that the formalism applies to arbitrary curvilinear coordinates. 
Various schemes exist to minimize the residual function R{x)] in the Galcrkin 
method one requires that there be no overlap between the residual and an 
arbitrary B-spline function 



v{x)dxBk{x)R{x) = 



(21) 



This so called Galerkin condition amounts to a global reduction of the residual. 
Applying the Galerkin condition and inserting the B-Spline expansions for 
f{x) and g{x) results in 



E 



v{x)dxBkix)OBi{x) 



E 



v{x)dxBk{x)Bi{x) 



b' = . (22) 



Defining the matrix elements 



Oki = j v{x)dxBk{x)OBi{x) , Gki = j v{x)dxBk{x)Bi{x) (23) 
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transforms the (differential) operator equation into a matrix x vector equation 



which can be implemented on modern vector or parallel computers with high 
efficiency. The matrix Gki is sometimes referred to as the Gram matrix; 
it represents the nonvanishing overlap integrals between different B-Spline 
functions (see Fig. ||). We eliminate the expansion coefRcients a', 6' in the 
last equation by introducing the function values at the lattice support points 
Xa including both interior and boundary points. 

The upper (U) and lower (L) components of the spinor wavefunctions 
defined earlier are represented on the 2-D lattice (r^jZ^) by a product of 
Basis Splines Bi{x) evaluated at the lattice support points 



We are also extending our previous B-spline work to include nonlinear grids. 
Use of a nonlinear lattice should be most useful for loosely bound systems 
near the proton or neutron drip lines. Non-Cartesian coordinates necessitate 
the use of fixed endpoint boundary conditions; much^ort has been directed 
toward improving the treatment of these boundaries Ej. 

4 Numerical tests and results 

We expect our Spline techniques to be superior to the traditional harmonic 
oscillator basis expansion method in cases of very strong nuclear deforma- 
tion. To illustrate this point, we have performed a numerical test using a 
phenomenological (Woods-Saxon) deformed shell model potential. We calcu- 
late the single-particle energy spectrum for neutrons in *'^Ca for quadrupole 
deformations ranging from strong oblate {(32 = — f.25) to extreme prolate 
(/32 = -1-2.25). The results are shown in Fig. ^. Apparently, for /32 = 
we correctly reproduce the spherical shell structure of magic nuclei. As (32 
approaches large positive values our s.p. potential approaches the struc- 
ture of two separated potential wells; as expected, we observe pairs of levels 
with opposite parity that are becoming degenerate in energy. The largest 
quadrupole deformation corresponds physically to a symmetric fission config- 
uration. Clearly, such configurations cannot be described in a single oscillator 
basis, which confirms the numerical superiority of the B-Spline lattice tech- 
nique. 
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Figure 3. Single-particle neutron levels for IgCa as function of quadrupole deformation in 
the Woods-Saxon shell model. Solid lines indicate positive parity and broken lines negative 
parity levels. 



In a second test calculation, we have investigated the properties of a nu- 
cleus near the neutron drip line. During the last decade the discovery of a 
'neutron halo' in several neutron-rich isotopes generated a great deal of in- 
terest in the area of weakly bound quantum systems. The halo effect was 
first observed in ^|Li, which consists of three protons and six neutrons in a 
central core and two planetary neutrons which comprise the halo. By adjust- 
ing the depth of the Woods-Saxon potential so that the separation energy of 
the last bound neutron is only 10 keV, i.e. very close to the continuum, we 
were able to determine this neutron wavefunction on the lattice which shows 
a very large spatial extent (see Fig. We conclude that the B-Spline lat- 
tice techniques are well-suited for representing weakly bound states near the 
drip lines; a similar calculation in the basis expansion method would require 
a large number of oscillator shells. 

We now discuss our numerical results for the selfconsistent Hartree-Fock 
calculations with Skyrme-M* interaction and BCS pairing. This is a special 
case of the HFB equation with a constant pairing matrix element. In Fig. 
H we display the proton density for a heavy nucleus, ^||Gd, calculated with 
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Figure 4. (a) Woods Saxon potential for ^^^Li; (b) spin- up component of the wave function 
for the last bound neutron. 




BCS pairing code; (b) measured charge distribution. 



our new 2-D (HF+BCS) code. It should be noted that ALL 154 nucleons are 
treated dynamicaUy (no inert core approximation). The theoretical charge 
density looks quite similar to the experimental result which is shown on the 
right hand side. 

For several spherical nuclei, we have also compared the selfconsistent s.p. 
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Table 1. Total binding energy and s.p. energy levels for O, calculated in the Hartree- 
Fock + BCS pairingi-feteory with the SkM* force. Comparison of results from ID radial 
finite difference code |£0 and our new 2D Spline-Galerkin code. We obtain the same level of 
accuracy despite the 25 times larger lattice spacing A = 0.625 fm. 



ID Radial 2D Spline-Galerkin 
A = 0.025fm A = 0.625fm 



Etot -127.73 MeV -127.48 MeV 

Esi/2in) -33.31 MeV -33.29 MeV 

-Ep3/2(n) -19.88 MeV -19.86 MeV 

-Epi/2W -13.55 MeV -13.53 MeV 

-Esi/2(p) -29.74 MeV -29.72 MeV 

Ep3/2ip) -16.48 MeV -16.45 MeV 

"pi/2(p) -10.27 MeV -10.26 MeV 



E. 



energy levels of our 2-D Spline-Galerkin code with a fully converged 1-D radial 
calculation. The result is shown in Table 1. 

Plans and Future Directions 

Having validated our new (HF-I-BCS) code on a 2D lattice with the Spline- 
Galerkin method, we plan to proceed as follows: We are currently working 
on the 2D HFB implementation with a pairing delta-force. After that, we 
will generalize the code utilizing the full SkP force with mean pairing field 
and pairing spin-orbit term. We will also add appropriate constraints, e.g. 
Q20: Qsoj i^Ja; for Calculating potential energy surfaces and rotational bands. 
As we compare the observables (e.g. total binding energy, charge radii, densi- 
ties pp „(r), separation energies for neutrons and protons, pairing gaps) with 
experimental data from the RIB facilities, it will almost certainly be necessary 
to develop new effective N-N interactions as we move farther away from the 
stability line towards the p/n drip lines. 
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